package scu.maqiang.elasticity;

import scu.maqiang.fes.BVPType;
import scu.maqiang.fes.FES3T42;
import scu.maqiang.mesh.Mesh3T4;
import scu.maqiang.mesh.MeshFormat;
import scu.maqiang.mesh.ScalarFunc;
import scu.maqiang.numeric.Direct;
import scu.maqiang.numeric.MVO;
import scu.maqiang.numeric.NewIterSSolver;
import scu.maqiang.numeric.SRMatrix;

public class ElasticityPressure3DT4 {

	public static void main(String[] args) {
		String meshFile = "MeshFile\\NewCantilevel3DT.msh";
		Mesh3T4 mesh3D = new Mesh3T4().readMesh(meshFile, MeshFormat.msh);
		int dispLabel = 10;
		ScalarFunc dispRegion = (xyz, label, param) -> Math.abs(xyz[0]) < 0.00001? 1.0:-1.0;
		mesh3D.setBoundaryLabel(dispRegion, null, dispLabel);
		int forceLabel = 20;
		ScalarFunc forceRegion = (xyz, label, param) -> Math.abs(xyz[0]) > 79.99999? 1.0:-1.0;
		mesh3D.setBoundaryLabel(forceRegion, null, forceLabel);
		mesh3D.toTecplot("cantilevel.dat");
		
		double E = 1.0e6, Nu = 0.3;
		FES3T42 fs = new FES3T42(mesh3D);
		int nv = mesh3D.getNv();
		int dofElastic=  fs.GetNdof();
		SRMatrix K = new SRMatrix(dofElastic);
		double[] RHS = new double[dofElastic];
		fs.assembleStiff(new double[]{E, Nu}, BVPType.COMMON, K);
		fs.assemblePressure(new double[] {-100.0}, BVPType.COMMON, RHS, forceLabel);
		fs.applyBC_MBN(K, Direct.All, dispLabel);
		double[] x = new double[dofElastic];
		fs.applyBC_MBN(RHS, Direct.All, 0.0, dispLabel);
		NewIterSSolver solver = new NewIterSSolver(K);
		solver.setEps(1.0e-16);
		solver.PCGSSOR(RHS, x, 1.5, 2);
		double[][] disp = new double[3][nv];
		fs.extractComponent(x, disp);
		System.out.println("Displacement: ");
		System.out.println("u min: " + MVO.min(disp[0])+ " \tmax: " + MVO.max(disp[0]));
		System.out.println("v min: " + MVO.min(disp[1])+ " \tmax: " + MVO.max(disp[1]));
		System.out.println("w min: " + MVO.min(disp[2])+ " \tmax: " + MVO.max(disp[2]));
		
		double[][] strain = new double[6][nv];
		double[][] stress = new double[6][nv];
		double[] eqStress = new double[nv];
		fs.computeStrainStressEqStress(BVPType.COMMON, new double[] {E, Nu},  x, strain, stress, eqStress);
		System.out.println("Strain: ");
		System.out.println("ex min: " + MVO.min(strain[0])+ " \tmax: " + MVO.max(strain[0]));
		System.out.println("ey min: " + MVO.min(strain[1])+ " \tmax: " + MVO.max(strain[1]));
		System.out.println("ez min: " + MVO.min(strain[2])+ " \tmax: " + MVO.max(strain[2]));
		System.out.println("exy min: " + MVO.min(strain[3])+ " \tmax: " + MVO.max(strain[3]));
		System.out.println("eyx min: " + MVO.min(strain[4])+ " \tmax: " + MVO.max(strain[4]));
		System.out.println("exz min: " + MVO.min(strain[5])+ " \tmax: " + MVO.max(strain[5]));
		
		System.out.println("Stress: ");
		System.out.println("tx min: " + MVO.min(stress[0])+ " \tmax: " + MVO.max(stress[0]));
		System.out.println("ty min: " + MVO.min(stress[1])+ " \tmax: " + MVO.max(stress[1]));
		System.out.println("tz min: " + MVO.min(stress[2])+ " \tmax: " + MVO.max(stress[2]));
		System.out.println("txy min: " + MVO.min(stress[3])+ " \tmax: " + MVO.max(stress[3]));
		System.out.println("tyx min: " + MVO.min(stress[4])+ " \tmax: " + MVO.max(stress[4]));
		System.out.println("txz min: " + MVO.min(stress[5])+ " \tmax: " + MVO.max(stress[5]));
		
		System.out.println("Von Mises Stress: ");
		System.out.println("min: " + MVO.min(eqStress)+ " \tmax: " + MVO.max(eqStress));
	}

}
